%% Script for synthesizing CGXHs for complex-amplitude control and simulation of FTH experiments 
% Double-Constraint Gerchberg-Saxton CGH Synthesis of a CGXH for the visible light experiments;
% 
% Written by Kahraman Keskinbora kahraman@mit.edu keskinbora@is.mpg.de
% in 2021 for the paper "Maskless Fourier Transform Holography"
%
% References:
% Chang et al., Appl. Opt. 54, 6994-7001 (2015)
% Voelz, Computational fourier optics: a MATLAB tutorial (SPIE, 2011)
%
%
close all;
clear all;
clc;

load cmap.mat

%% Units;

meters      =1;
microns     =1e-6 * meters;
nanometers  =1e-9 * meters;
cm          =10^-2 * meters;

seconds     =1;
miliseconds =1e-6 * seconds;
nanoseconds =1e-9 * seconds;

hertz       =1/seconds;
kilohertz   =1e6 * hertz;
gigahertz   =1e9 * hertz;
terahertz   =1e12* hertz;
petahertz   =1e15* hertz;

eV          =1;
keV         =eV * 1000;
meV         =eV / 1000;


%% Constants;

c0          = 299792458 * meters/seconds;
e0          = 8.854187817e-12 * 1/meters;
u0          = 1.2566370614e-6 * 1/meters;

%% Target Definition

L           = 2*cm;          % Initial Sim bounding box
dx          = 8*microns;          % Step size
Nx          = L/dx;
x           = -L/2:dx:L/2-dx;
y           = x;
[X,Y]       = meshgrid(x,y);
lambda      = 635*nanometers;
f           = 30*cm; 
seperX      = -500*microns;
seperY      = -500*microns;

% Rectangular Sample Beam
fov_FW      = 500*microns;                    % Sample beam FOV width               
fov_hw      = fov_FW/2;                     % FOV half-width
Fov_spot    = rect((X+seperX)./(2*fov_hw)).*rect((Y-seperY)./(2*fov_hw));

% Reference Beam
Ref_spot_D  = 8*microns;
Ref_spot    = circle(Ref_spot_D,0,0,x);

IR          = 150;                          % Intensity ratio of reference spot to sample beam
target_abs  = Fov_spot+IR.*Ref_spot;
target_phase= (Fov_spot+Ref_spot);
target      = target_abs.*target_phase;

CGXH_w      = 0.9660*cm;                      % CGXH size
support     = rect(X./CGXH_w).*rect(Y./CGXH_w);     % CGXH support

%% Double-Constraint Gerchberg-Saxton Algorithm
% Initate the Input Field;

G=zeros(Nx);

PWI         = 1;                                  % Plane wave illumination

iterations  = 20;
source      = G;
counter     = 1;
tc          = 0;

% DC-GSA

alpha   = 1;
beta    = 0.01;
gamma   = 0.1;
Npix    = Nx.^2;  
mask1   = abs(target);
mask2   = (1-abs(target));
z       = f;


 for i=1:1:iterations           
    tic;
    A=support.*exp(1i*angle(source));
    B=propRS2(A,L,lambda,z);
    C1=(alpha.*abs(target)-beta.*abs(B)).*exp(1i.*angle(target));
    C2=(gamma.*abs(B)).*exp(1i.*angle(B));
    C=C1+C2;
    D=propRS2(C,L,lambda,-z);
    source=D;
    t(i)=toc;
    tc=tc+t(i);

    RMS(i)=error_function_Miao(abs(B),abs(target))
    figure(1)
    plot((RMS))
    drawnow

 end

RMS(i)

%% Export the synthesized hologram for displaying on the SLM screen
holo        = angle(A)+pi;

holo2       = holo./max(max(holo));
DC          = 0.5;
binary_CGH  = ((0).*(holo2<DC)+1.*(holo2>=DC)).*support;

% Crop to dimensions
binary2     = binary_CGH(649:1848,649:1848);

% Export bitmap to be displayed on the SLM screen
imwrite(binary2,'SLM_StIXH_TEST_bin.png')
